In granular assemblies, there are inhomogeneous contact forces and geometrical structures on the scale of granular particles, which makes the granular system show complicated properties such as nonlinearity, disorder, and anisotropy,^[1] and the macroscopic description of the material behavior is very difficult. The inhomogeneous stress distribution in quasistatic granular packings will cause the contact force between some particles to be greater than the average stress; these particles usually exist in the form of chains, and support the whole system.^[2] It is still not clear that how the heterogeneity of the stress distribution and the existence of the force chain affect the properties of the granular system, whether or not the granular system can be regarded as an elastomer. If it can be regarded as an elastomer, it is also controversial whether the elastic theory is applicable to describe the granular system at various scales.^[2] To obtain the microstructure evolution and macroscopic mechanical properties of a system, the most intuitive research methods are the photoelastic experiment,^[3] x-ray,^[4] and nuclear magnetic imaging,^[5] which are expensive. Moreover, these methods also have some difficulties to detect the real granular system. Besides these methods, acoustic propagation experiments can be used to understand the relationships among the micro-stress distributions, geometrical structures, and macroscopic properties. The nonlinearity and disorder of the granular system play important roles in the propagation of mechanical energy, and the microstructure and stress distribution correspond to different propagation modes.^[6] The acoustic propagation experiment provides an effective way to detect the internal structure of the real granular system. For example, the low frequency acoustic detection experiment can determine the nonlinear elastic property and structural anisotropy of the granular media,^[7–9] while short-wave multiply scattered sound waves enable one to detect tiny changes of the contact network configuration at the grain level,^[10,11] as well as the use of acoustic detection techniques to study seismic triggering^[12] and the formation of shear bands, etc.^[13]

In the acoustic experiment, the determination of the sound velocity is usually done by the time of flight method.^[14,15] This method obtains the time and the corresponding velocity of the sound wave in the granular system through the real path because sound waves favor speed and amplitude along force chains.^[16,17] The waveform and velocity obtained by the time of flight method reflect information such as stress, strain, stiffness tensor, and attenuation along the wave propagation path. There is another method to obtain the sound velocity apart from the time of flight method. If the granular system is regarded as a continuous medium, the elastic potential energy of the system can be calculated according to the contact relation based on the effective medium theory (EMT)^[8] or granular solid hydrodynamics (GSH) theory.^[15] Then, the stress can be obtained by taking partial derivative of the potential energy with respect to strain, and the sound velocity can be obtained from the fourth-order stiffness tensor.^[18,19] EMT and GSH theories usually deal with homogeneous cases without considering the inhomogeneity of geometry and stress distribution, and the given stiffness tensor and sound velocity are the average effect of the whole system.^[15]

Whether from the time of flight method or the stiffness tensor, the sound velocity is the response of the system stress structure. For example, the relationship between the normal stress p and sound velocity of the disordered packing granular system obtained by the time of flight method v_tof meets v_tof ∝ p^1/6 while this dependence is at higher pressures, but the behavior at lower confining pressures is more like v_tof ∝ p^1/4.^[7,15,20] However, the sound velocity calculated from the stiffness tensor v_c meets the relation v_c ∝ p^1/6.^[2,8,21,22] The compression waves are more sensitive to the stress-induced anisotropy, whereas the shear waves are more sensitive to the fabric anisotropy.^[9] The two-dimensional photoelastic experiment shows that the sound velocity v has a nonlinear behavior with the applied force F as v = KF^α, where the exponent α depends on the contact number of the grains, while the prefactor K brings information about the anisotropy of the medium.^[17,23] Therefore, the study of the relationship between sound velocity and the internal structure of particles can not only detect a change in the structure of the granular system, but also recognize its nature, and it is of great interest when trying to predict abrupt events, namely, avalanches, earthquake, or unjamming event.^[16,24]

The anisotropic granular system^[1] is obtained experimentally with different preparation methods^[9] (preparation protocols of rain deposition, decompaction) or different stress loading methods; however, it is difficult to measure the internal stress distribution of the system. Usually, the force on the boundary wall is used to substitute the stress distribution in the system,^[9] and the calculation of stiffness tensor has also been carried out by some simplification and approximation.^[14] The limitation of experimental means makes it difficult to accurately describe the particle system. The acoustic experiment, especially the acquisition of the stable shear wave, is also more difficult and higher requirement for equipment is needed.

The discrete element method (DEM) can obtain the force and deformation of each particle, and simulate the sound wave through the simple harmonic vibration or simple harmonic force of the boundary wall, making it convenient to obtain the waveform and sound velocity. Furthermore, DEM can obtain the required physical quantities of EMT and GSH theories to explain acoustical properties of granular systems, for example, the coordination number used to explain the relation v ∝ p^1/4,^[25] the volume and shear modulus used to calculate the sound velocity,^[8] and so on.

In this paper, four different packing structures (S45, H60, S90, and D) in the two-dimensional photoelastic experiment are simulated by DEM.^[17] The relationship between compressional and shear wave velocities and normal stress is calculated by the time of flight method. The influence of different packing structures on the sound velocity under the same normal stress is analyzed. The stress tensor and the stiffness tensor under different packing structures are calculated. The wave velocities obtained by the stiffness tensor are compared with those by the time of flight method. The variation of sound velocity with stress and packing structures can be explained by the difference of system stiffness and stress tensor. The results are consistent with the existing experiments and theories.^[17] The change of sound velocity reflects the difference of stress distribution and geometrical structure in the system, which can provide a reference for acoustically probing the structure of system.

Four kinds of packing structures used to calculate the sound velocity are shown as S45, H60, S90, and D (disordered) in Fig. 1. The particle numbers of the systems, in turn, are 4858, 4865, 4900, and 4165.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) Different packing structures.

The contact between particles is simulated by the simplified Hertz–Mindlin–Deresiewicz model. The force and torque expressions are listed in Table 1.

Table 1.

Table 1.

Table 1. Force and torque between particles. . Property Equation Normal spring force Normal viscous force Tangential spring force Tangential viscous force Torque qi = ric × (fct + fdt) Effective Young’s modulus E* = E/(2 − 2v2) Tangential displacement vector Tangential maximum displacement δt, max = μs(2 − v)/(2 − 2v)δn Relative sliding velocity vij = vi − vj + ωi × ric − ωj × rjc Normal relative velocity Tangential relative velocity vt,ij = vij − vn,ij

Table 1. Force and torque between particles. .

The parameters of the simulations are listed in Table 2.

Table 2.

Table 2.

Table 2. Parameters of the simulations. . Parameters Value SI-units Length-unit 1 1 mm Mass-unit 1 1 mg Time-unit 1 1 ms Size of system 70 × 70 70 mm × 70 mm Particle density 2.65 2.65 ×103 kg/m3 Particle mass 2.08 2.08 mg Moment of inertia 0.26 2.6 × 10−11 kg·m2 Young’s modulus 107 10 GPa Background damping 700.0 0.7 kg/s Friction coefficient 0.40 0.40 Poisson’s ratio 0.20 0.20 Normal viscosity coefficient 0.20 0.20 Tangential viscosity coefficient 0.20 0.20 integration time-step 2.0 ×10−4 2.0 ×10−7 s

Table 2. Parameters of the simulations. .

Using walls as boundary conditions, the top wall is subjected to normal stress, while the other three walls are fixed in the loading process. When the system achieves equilibrium completely under a certain normal stress, the sound wave is simulated by a simple harmonic vibration of the bottom wall. The simple harmonic vibrations in the vertical and horizontal directions represent the compression and shear waves, respectively, as shown in Fig. 2. The vertical and horizontal vibration frequencies are both 20 kHz. The waveform is determined by comparing the force of the top wall when the sound wave exists with the force without sound wave. Figure 3 shows the input and output waveforms of the disordered structure under a normal stress of 3 MPa.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Sound wave model.

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) Input and output waveforms: (a) input waveform, (b) output compression waveform, and (c) output shear waveform.

Taking the time at 1/10 of the first wave peak as the flight time of the sound wave, the distance between the top and bottom walls is divided by the time to give the velocity. As shown in Fig. 3, the peak value of the compression wave is 6.7377 × 10⁻³, and 1/10 of the value corresponds to the time t = 9.75 × 10⁻⁵ s. The distance between the top and bottom walls is l = 60.455 mm, and the velocity is v_tof = 620.05 m/s, where the subscript tof indicates the velocity obtained by the time of flight method.

The velocities by the time of flight method and stiffness tensor are related to the packing structures. The compression wave velocities of H60 and S90 are larger than those of S45 and disordered structure under the same normal stress, and the relative velocity difference can reach 30%. For instance, at the normal stress of 3 MPa, m/s is nearly the same as , bigger than m/s, and much bigger than m/s. The difference between compression wave velocities of different structures reveals that the compression wave velocity is sensitive to the volume fraction. The volume fraction of H60 (hexagonal closest packing) is bigger than that of D, and the volume fraction of S45 is the smallest, which means that a bigger volume fraction corresponds to a bigger compression wave velocity. Ignoring boundary conditions, the volume fraction of S90 is the same as that of S45, but the compression wave velocity of S90 is nearly the same as that of H60 and bigger than that of S45, which is due to the special structure of S90 only for the propagation of compression wave as illuminated by the following discussion.

The relationship between shear wave velocities and packing structures are different. The shear wave velocity of S45 is nearly the same as that of H60, and bigger than that of disordered packing. For example, at the normal stress of 50 kPa, m/s is bigger than m/s. For the structure of S45, the compression wave velocity of which is the smallest, the shear wave velocity is bigger than that of disordered packing. The reason is that the shear and compression wave velocities of S45 are the same because of structural symmetry, while the shear wave velocities of D and H60 are smaller than their compression wave velocities. This phenomenon reveals that structure has a greater impact on the shear wave velocity than the volume fraction. In other words, shear wave velocity is more sensitive to the packing structure, which is consistent with the results in Refs. [9] and [17].

The velocity difference between various structures can also be explained by different stiffness tensors. Different structures mean different geometries and force distributions, and then determine the stiffness tensors. From the ideal of a continuum, velocity depends on stiffness, and in fact, the velocity difference is the stiffness difference determined by the geometry. The change in stiffness determines the change in velocity, especially for greater normal stress (more than 1 MPa).

In addition to stiffness, the difference of velocity can be interpreted by the dispersion of elastic waves in granular phononic crystals. In the case of hexagonal close packed structure, as shown in Fig. 6, the contacts are modeled by the normal, tangential, and rolling springs with the stiffnesses k_n, k_t, and k_r. The position of the particle in the center of the lattice (m,n) is (md, nd), where d is the diameter of the particle, and positions of other particles are shown in Fig. 6(a). The kinetic behavior of the phononic lattice can be deduced from the Lagrangian L of the particle (m, n), which reads

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Contact model of hexagonal close packed structure: (a) normal spring, (b) tangential spring, and (c) rolling spring.

The dispersion relation gives the velocity of the shear wave and that of the compression wave, in the long-wave limit:

For the square phononic crystal (S90), the compression wave velocity is^[30]

The ratio of compression wave velocities for H60 and S90 is , and the result is coinciding with Figs. 4(b) and 4(c). The ratio of the compression wave velocity and the shear wave velocity for H60 is , which is bigger than the ratio (≈ 1.3) calculated by Eq. (3).

For S45, because of the symmetry, the first two equations of motions are the same, and the velocities of the compression and shear waves should be equal, which is in agreement with Figs. 4(a) and 5(a).

The velocities from the lattice dynamics are ideally given, and the system is always symmetric and homogeneous. However, for a realistic granular system, it is difficult to prepare a sample in which the force distribution is the same as its geometry even by simulation. That is the reason for the deviation between the velocities by the time of flight method and those by the stiffness or lattice dynamics.

The relationship between the normal stress and velocity derived from the stiffness tensor is v_c ∝ p^1/6 as shown in Figs. 4 and 5, which is determined by the Hertz model and is consistent with the existing theory.^{[2,8,15,21,22]} When the normal stress p is larger than 1 MPa, there is minor difference between the sound velocity obtained from the time of flight method and those from the stiffness tensors. When the normal stress p is less than 1 MPa, the velocity from the time of flight method is less than that from the stiffness tensor, that is, when the velocity is under a smaller normal stress, the velocity meets v_tof ∝ p^α, where α > 1/6, and the linear correlation coefficient is greater than 0.95 in the fitting simulation result with α = 1/4 obtained by experiments.^[7,20]

The shear wave velocities by the time of flight method have an overall deviation from the velocities (Fig. 5(a)) by the stiffness tensor, which may be related to the internal stress distribution. Figures 7(a) and 7(c) show the horizontal and vertical force distributions in the sample (sample 1) which is used to calculate Fig. 5(a) with the normal force p = 400 kPa. The force distribution of sample 1 is so inhomogeneous that the horizontal force has the layer structure. The difference between the average stiffness tensor and the stiffness tensor in the actual propagation path of the acoustic wave is relatively large, resulting in difference of the sound velocities obtained by these two methods.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Force distributions in samples of different preparation methods: (a) Fx in sample 1, (b) Fx in sample 2, (c) Fy in sample 1, (d) Fy in sample 2.

Figures 7(b) and 7(d) show the force distributions in sample 2. Because of the small volume fraction of S45, particles are more prone to deviate from the setting geometry, and sample 2 is prepared by increasing the damping coefficients c_n and c_t to balance the deviation. The shear velocities of different damping coefficients are shown in Fig. 8. It can be seen that the pattern of shear velocity tends to keep accordance with the patterns of H60 and D, and the transformation from v_tof ∝ p^1/4 to v_tof ∝ p^1/6 still exists. In Ref. [17], the transformation is explained by the fact that the force-deformation curve of the two contacts geometry shows a different power law from that of the four contacts geometry. From the point of the coordination number, Reference^[25] also explains the reason of the change. The Hertz contact model used in this work keeps the relationship f ∝ δ^3/2 unchanged, and the coordination number remains unchanged for S45, H60, and S90; in this case there is still v_tof ∝ p^1/4 which means the cause is more complicated. An explanation for the transformation is the dissipation caused by the disturbance of the wave propagation as suggested by Ref. [15].

Fig. 8.

	Figure Option ViewDownloadNew Window
	Fig. 8. (color online) Influence of damping on the shear velocity: (a) Cn = 0.75, Ct = 0.5, and (b) Cn = 0.5, Ct = 0.25.

The 1/6 exponent of Hertz scaling is limited to homogeneous and non-dissipative samples, however, dissipation is the nature of reality in the particle system. When the particle system achieves balance under external force, the strains and the relative positions between particles are unchanged. In the state of equilibrium, the system can be considered as an example of elastomer. When the system is disturbed by waves, the relative velocities and positions between the particles emerge. As a result of the relative velocity and position, the force chains in the system have the possibility of relaxation, then the dissipation emerges, which causes the system to deviate from elastomer, that is, deviate from the 1/6 exponent of Hertz scaling.

The fluctuation of energy under different normal stresses during the shear wave propagation is shown in Fig. 8. The elastic potential energy of the system is easy to obtain by integrating f_cn and f_ct against δ_n and δ_t, respectively, for all contacts, because the relative velocities between particles are all zero, and the potential energy is the total energy of the system. The energies of the input waves of different stresses are equal to 2.4 × 10⁻⁷ J, and the kinetic energy of particles, which is zero before propagation, changes as Fig. 9(a), with the ratio of kinetic energy to elastic potential energy shown in Fig. 9(b).

Fig. 9.

	Figure Option ViewDownloadNew Window
	Fig. 9. (color online) Fluctuation of energy during wave propagation: (a) kinetic energy of particles, and (b) ratio of kinetic and potential energy.

For small normal stress p = 100 kPa, the kinetic energy of particles is larger than the input wave energy, which means the difference of the energy comes from the initial elastic potential energy of the system. That is to say, some of the elastic potential energy stored in force chains is dissipated by the wave disturbance, especially in the tangential direction, and the number of the broken tangential springs at every time step is shown in Fig. 10. The strength of force chains is different from the initial elastic state, which probably cause the deviation from the 1/6 exponent law.

Fig. 10.

	Figure Option ViewDownloadNew Window
	Fig. 10. (color online) The number of tangential sliding, p = 100 kPa.

In Fig. 9(b), with the normal stress increasing, the kinetic energy of particles during the wave propagation gets smaller, especially for stresses larger than 1 MPa; the ratio of kinetic energy to potential energy is very close to zero, that is, the wave disturbance almost has no influence on the initial elastic state, and figure 9(b) echoes the 1/6 exponent of Hertz scaling for p > 1 MPa.

From another point of view, in the contact model used in this paper, c_n and c_t are used to balance the relative velocity and position between the particles. The larger the values of c_n and c_t, the more difficultly the particle moves; in other words, for the static system, c_n and c_t contribute to prevent the dissipation caused by the relative motions of the particles, that is why the shear wave velocity tends to keep accordance with the 1/6 exponent of Hertz scaling as c_n and c_t mounted, as shown in Figs. 5(a) and 8. For large normal stress or big volume fraction, the relative motions of particles are difficult, and the system is also difficult to deviate from elastomer, so the 1/6 exponent of Hertz scaling is easy to satisfy.

In the end, velocities of four packing structures reflect different geometries and force distributions; the compression wave velocity is sensitive to the volume fraction of the structure, and the shear wave velocity is more sensitive to the geometrical structure itself. For the stress dependence of the sound wave velocity, the dissipation is the probable cause of the deviation from 1/6 exponent of Hertz scaling, and further research will be done in the future.